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ABSTRACT 

We follow the sinking of two massive black holes in a spherical stellar system 
where the black holes become bound under the influence of dynamical friction. 
Once bound, the binary hardens by three-body encounters with surrounding 
stars. We find that the binary wanders inside the core, providing an enhanced 
supply of reaction partners for the hardening. The binary evolves into a highly 
eccentric orbit leading to coalescence well beyond a Hubble time. These are the 
first results from a hybrid "self consistent field" (SCF) and direct Aarseth N- 
body integrator (NBODY6), which combines the advantages of the direct force 
calculation with the efficiency of the field method. The code is designed for 
use on parallel architectures and is therefore applicable to collisional iV-body 
integrations with extraordinarily large particle numbers (> 10 5 ). This creates 
the possibility of simulating the dynamics of both globular clusters with realistic 
collisional relaxation and stellar systems surrounding supermassive black holes 
in galactic nuclei. 



1. Introduction 

Currently the standard picture of galaxy formation involves the collapse of baryonic 
matter in hierarchically clustering dark matter halos and the subsequent building of big 
galaxies from small ones via merging processes e.g., (Peebles 1993; Diaferio et al. 1999; 
Kauffmann et al. 1999a,b). While recent cosmological simulations can adequately reproduce 
many global properties of galaxies and their correlations, the details are still very much 
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dependent on the gas physics and stellar feedback involved (see e.g., Navarro and Steinmetz 
(2000)). Additionally, most, if not all, galaxies harbor supermassive black holes in their 
center (Magorrian et al. 1998; Richstone et al. 1998; Kormendy and Richstone 1995). 
Correlations have been recently detected between black hole masses, galaxy masses, and 
central velocity dispersions in galaxies (Ferrarese and Merritt 2000; Gebhardt et al. 2000). 
These correlations are strong evidence that black holes in galactic nuclei are linked to the 
dynamical history of their host galaxies. Haehnelt and Kauffmann (2000) and Kauffmann 
and Haehnelt (2000) demonstrate how this is consistent with the framework of semi-analytic 
models that follow the formation and evolution of galaxies in a cold dark matter-dominated 
universe. They assume supermassive black holes are formed and fueled during major mergers, 
qualitatively explaining many aspects of the observed evolution of galaxies, including the 
observed relation between bulge luminosity, velocity dispersion, and central black hole mass. 
As already discussed by Begelman et al. (1980), such a scenario requires the formation 
of galactic nuclei containing at least two black holes, depending on the black hole merger 
rate relative to the galaxy merger rate. However, there is very little observational evidence 
for massive black hole binaries (Lehto and Valtonen 1996; Halpern and Eracleous 2000). 
This conflict between theory and observations has become known as the "sinking black hole 
problem" . As an alternative to minimally impacting stellar dynamical processes, Gould and 
Rix (2000) and Armitage and Natarajan (2002) have proposed mechanisms which lead 
to rapid decay of massive black hole orbits and subsequent black hole mergers in galactic 
centers. Also, Begelman et al. (1980) offered the solution that gas accretion could dominate 
the orbital decay in the intermediate phase of the sinking black hole problem when dynamical 
friction becomes inefficient. However, as we will discuss later, dynamical friction, as laid out 
by Chandrasekhar (1943), is not sufficiently effective by itself to lead to rapid coalescence 
of black hole binaries. 

If there are no quick mergers, multiple black hole nuclei could lose black holes through 
slingshot ejections (Valtonen et al. 1994). Once a binary system becomes hard, the high 
orbital velocities of the black holes allow further hardening through close encounters and 
three-body interactions with stars. Such processes will evacuate field stars from the sur- 
roundings of the binary, therefore it can be argued that the stellar scatterings cannot produce 
rapid coalescence. The preceding argument assumes that the center of mass of the binary 
does not move with respect to the stellar system. However, we will show that even with a 
fairly symmetrical initial setup the binary gains some linear momentum. This introduces a 
wandering motion which exceeds the expectations from equipartition. The wandering of the 
binary guarantees an adequate supply of stars for binary hardening and rapid coalescence 
through purely stellar dynamical processes. 

Our new computational method allows us to study in detail three-body interactions of a 
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black hole binary with field stars. Although one may argue that the perturbing mass of the 
field stars is small compared to the black hole mass and should have negligible impact, there 
are many stars, and each encounter can lead to changes in binding energy and eccentricity 
of the black hole binary. In fact, our models show that the black hole binary keeps a rather 
high eccentricity due to the encounters. Thus high eccentricity will speed up gravitational 
radiation mergers very efficiently, and is, as noted by Gould and Rix (2000) and Armitage 
and Natarajan (2002), a way to expedite massive black hole mergers in a purely stellar 
dynamical way. 

The correct theoretical prediction of the frequency of black hole mergers in galactic 
environments will be important in the search for gravitational waves. The merging of super- 
massive black holes of 3 x 10 4 to 3 x 1O 7 M in the nuclei of merging galaxies and protogalaxies 
can be detected with high signal-to-noise at redshifts from < z < 100 (Phinney 2000) by 
the Laser Interferometer Space Antenna (LISA) (Danzmann 2000). 

Previous attempts to quantify this prediction have been made by either solving the 
perturbed two and three-body problem in simplified models (Mikkola and Valtonen 1992), 
direct iV-body models (Makino et al. 1993; Makino 1997), or a combination of the two 
(Merritt and Quinlan 1998; Quinlan and Hernquist 1997). Simulating binary black hole 
hardening is extremely challenging, algorithmically and computationally. Since the mass 
differences between the black holes and the stars is so large, high particle numbers are 
required in order to model the relaxation processes around the black holes accurately. The 
simulations have used softened particles on special purpose computers (Makino et al. 1993; 
Makino 1997) or a hierarchical hybrid code in which all forces involving the black hole 
particles are Keplerian (Merritt and Quinlan 1998; Quinlan and Hernquist 1997). These 
schemes used particle numbers in the order of 10 4 . 

In this paper, we describe a new hybrid field-particle code which treats all particles with 
orbits crossing the central regions of the system with a high precision direct iV-body method 
appropriate for collisional stellar dynamics. All other particles are integrated using a field 
method. In order to adapt both parts of the hybrid code to each other, the field method 
(approximating the potential exerted by a set of particles by a series expansion, referred to 
here as "SCF") had to be upgraded to a fourth order Hermite integrator. This integration 
also uses the time derivative of the potential, as in modern direct iV-body codes. 

In the following sections some details of the sinking black hole problem are introduced. 
Section 2 introduces the integration software used for the numerical experiments described 
in this paper. Section 3 is devoted to a comparison between the new collisional code with 
a well used workhorse simulator in this field called NBODY6 (Aarseth 1993, 1996, 1999), 
using its parallel implementation NBODY6++ (Spurzem and Baumgardt 2001; Spurzem 
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1999). The application of the code to the sinking binary black hole problem is reported in 
section 4. 



2. Collisional stellar dynamics with EuroStar 

A numerical simulation of the hardening phase (until the massive black holes start to 
radiate gravitational waves) must be able to accurately follow three-body encounters. For 
this reason, the Keplerian potential should not be softened in the denser parts of the system. 
Computationally, the central part of the system would best be treated using a collisional 
integrator. The code must be able to integrate encounters leading to large angle deflections 
in an efficient way, while requiring neither too much computing time nor introducing energy 
errors. The overall iV-body integration does not need to be symplectic, but should keep 
the energy error as low as possible. On the other hand, in a system showing a core halo 
structure, the bulk of the stars in the halo move under the influence of the mean field of the 
whole cluster. The halo part of the central galactic cluster can be integrated with a mean 
field method. 

In the new method (which we will refer to as EuroStar), both the collisional code 
NBODY6++ (Aarseth 1993, 1999; Spurzem and Baumgardt 2001) and the SCF method 
(Hernquist and Ostriker 1992; Hernquist et al. 1995; Zhao 1996; Sigurdsson et al. 1997; 
Holley-Bockelmann et al. 2001) are merged to optimize large- N collisional iV-body simula- 
tions. The star cluster, which is assumed to be in equilibrium, is divided into two sections 
by applying a critical angular momentum criterion. As shown in Figure 1, this allows for 
distinction between particles orbiting solely in the halo and ones which have trajectories 
leading through the core of the system. 

In a stationary gravitational point mass system, two-body relaxation leads to an ex- 
change between halo particles and core particles in a divided cluster. In a system of more 
than 10 4 particles, only a few particles cross the core halo border per dynamical time. This 
is very fortunate because it allows us to integrate the orbits of the halo particles with a 
collisionless method and the core particles with a collisional code. Exchanges of particles 
cause energy conservation problems since the contribution of a particle to the main potential 
would be changed from Keplerian to a sample point of a mean field in EuroStar. This means 
that switching of particles from the core to the halo part of the integrator and vice versa is 
not permitted. 

NBODY6+- 1- integrates trajectories of point masses in the core of the system. It is an 
Aarseth-type direct force integrator applying the Hermite integration scheme. NBODY6+- 1- 
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Fig. 1. — Schematic decomposition of a star cluster according to its angular momentum 
distribution. The hatched region symbolizes the distribution of angular momentum per unit 
mass as a function of radius. The solid curve above shows the angular momentum of particles 
on circular orbits with escape velocity. As such, the solid curve functions as an upper bound 
for the angular momentum at a given radius, or as a lower limit for the radius a star can 
reach with fixed angular momentum. If the cluster is divided into two parts by a critical 
value for /, as shown by the horizontal dashed line, halo particles which never reach the core 
can be distinguished from particles which pass deep into the core (vertical solid line). The 
dot-dot-dashed line shows the maximal angular momentum for a set of particles selected by 
an energy criterion. Such a selection would not affect all particles below that line. In fact, 
an energy criterion is not sufficient for selecting all collisional particles in our system. 
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gains its efficiency by implementing an Ahmad-Cohen neighbor scheme and individual block 
time steps (Ahmad and Cohen 1973; Aarseth 1999). Close interactions between particles 
are treated by regularization of the equations of motion (Kustaanheimo and Stiefel 1965). 
NBODY6++ scales well on parallel computer systems and can also be used with the GRAPE 
special purpose computer (Spurzem and Kugel 1999; Sugimoto et al. 1995). The Hermite 
scheme requires one to compute F, and F« at each time step where, 

f, = -£^. (1) 
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The relative distance between particles % and j is given by = r, — r,-. Accordingly, the 
relative velocity is = Vj — Vj. The extra effort of computing two direct force quantities 
allows one to approximate the particle's orbit to fourth order. By storing Fj and Fj from 
the previous time step it is possible to interpolate the next two higher derivatives and to 
apply a predictor /corrector scheme (Aarseth 1996). 

The SCF method qualifies for the collisionless part of a spherical system since then the 
basis functions are given analytically. This allows one to implement the Hermite scheme for 
SCF, which makes SCF an ideal far force extension to NBODY6++. Its drawback, however, 
is that this method restricts the input systems to have an approximately spherical particle 
distribution around the coordinate center. The possible asphericity depends on the number 
of spherical harmonics used for the potential expansion. In order to have better convergence, 
one-parameter basis functions for p„/ m (r) and $„; m (r) are used (Zhao 1996): 

Pmm(r) = V^—^^^—YU^ (3) 
r ( ~> (1 + r~) 2+a ( 2l+1 > 

®nlm(r) = ~ ~ FT ' Y lm (#, <p). (4) 
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The &n ] are called ultraspherical or Gegenbauer polynomials. The spherical harmonics are 
given by the Y^fi, ip). Once the A n i m for a certain set of particles are known, an analytic 
expression of the potential and the density is found. Due to the truncation of the expansion 
they represent the mean field and mean density. This means that the force at each position 
and its derivative can be computed using the following expressions, 

F(r) = -J2 A m m V<lw(r), (5) 

nlm 

^( r ) = -|(£ A>hnV<Wr))- (6) 

nlm 
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Since Equations (5) and (6) lead to a significant modification of the SCF scheme, we provide 
detailed form of these expressions in the Appendix. 



3. Testing the hybrid code 

As a first test for the new method, we have followed the last stages of an ongoing merger 
between two galaxies, each containing a supermassive black hole. For the initial setup, 
it is assumed that the stellar systems have already arranged themselves into a spherical 
system. The two, formerly central, supermassive black holes are moving through the stellar 
component with a speed on the order of the relative velocity between the two initial galaxies. 

In the present simulations, the stellar component is a realization of a Plummer model. 
The density and potential of the spherically symmetric Plummer model are given by (Plum- 
mer 1911): 

n 3M R 2 
P[r) ~ 4tt (i? 2 + r 2 )5/2' ^ 

$0) = -GM (8) 

^ (R 2 + r 2 y/ 2 ' 

The quantity M describes the total mass of the system and G is the gravitational constant. 
With the Plummer radius chosen to be R = 3n/16, the half mass radius of this system is at 
a radius of r h ~ 0.78 in the model units of the simulations. The total mass of the system 
M is set to unity. The gravitational constant G is set to unity as well, conforming to the 
model units described by Heggie & Mathieu (1986). The stellar system is centered around 
the origin. 

The black hole particles contain 1% of the system's total mass. Their initial positions 
are at x — ±0.5, and their initial velocities are 13.6% of the circular velocity at their initial 
radii. 

The black hole orbits are analyzed during the simulation assuming the orbit can be 
approximated by the classical two-body problem. The binding energies and eccentricities 
of the black hole orbits are computed from their relative distances and velocities, assuming 
a Keplerian potential. Once the black holes become bound, their two-body attraction is 
the most important force. The eccentricity of the binary e and the binding energy h are 
computed as follows, using the definition r = r a — r b . The vectors r a and r b denote the 
position vectors of the black holes a and b. From this it follows that (Boccaletti and Pucacco 
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1996), 
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This method of analysis provides a sensitive measure for the moment when the black holes 
become bound to one other. Furthermore, this way of analyzing the data also offers a precise 
tool for following the hardening of the binary. 

Since this sinking binary problem is the first application of EuroStar, its results are 
compared with those of the fully collisional code NBODY6++. For the comparison runs, 
16384 particles were simulated. Figure 2 shows the results for the two comparative runs. 
The plot on the left hand side in Figure 2 shows the eccentricity of the binary as a function 
of the simulated time in model units. The plot on the right hand side shows the two-body 
binding energy as a function of time. The binary becomes bound after 10 time units, in both 
cases. 

The fully collisional method and the hybrid code EuroStar show slightly different sinking 
rates for the binary at the beginning of the simulation. These differences result from the 
different density of collisional particles around the black holes in both codes. While the 
black holes in the fully collisional run suffer small angle encounters with every particle in 
the system, this is not possible in the hybrid code. Naturally, all particles treated by the 
mean field method can only interact with the system through changes in the mean potential. 
After the binary has become bound, the hardening process is driven by the stars which have 
a small enough impact parameter such that they have an encounter timescale smaller than 
the orbital timescale of the massive binary. Therefore, the hardening depends more on the 
two-body encounters with neighboring particles, which have large orbital velocities. 

With increasing simulation time, the binary locks into an oscillating motion around the 
center of mass of the stellar component. This motion does not extend beyond of the dense 
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Fig. 2. — Development of the orbital eccentricity of the black hole binary as a function of 
time in iV-body time units (= NBU, left graph) and its binding energy as a function of time 
(right). The results for the direct method are shown by the solid line (NBODY6++), the 
ones for the hybrid code (EuroStar) are given by the dashed line. Both methods use the 
same initial model with 16384 particles. Note, that if the binary is not yet bound, equations 
(10) and (11) formally yield values of a > and e > 1. This means that in that phase the 
black holes are still not yet gravitationally bound to one other. 
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galactic core. Effectively, the differences in the density of the collisional particles between 
the two methods vanish after the binary has become bound. Figure 2 reflects this by showing 
a parallel evolution of eccentricity and binding energy for times larger than 10 time units in 
the simulation. 



4. Hardening of a massive binary 

This new code is intended for raising the total particle number for collisional simulations 
of spherical iV-body systems. Hence the evolution of one, two or several massive bodies 
in a dense stellar cluster appears to be an ideal problem for EuroStar. This is why we 
are addressing here the problem of a sinking massive black hole binary in galactic centers. 
Another useful potential application os the dynamics of globular clusters. 

4.1. Initial conditions 

The particles representing the stellar component are distributed according to Plummer's 
model with R = 37r/16. The total mass of the stars is fixed at 0.98, while the black hole 
particles carry 0.01 each, so M = 1.0. This is a fairly high mass for the black holes compared 
to the total mass of the stellar system, since Ferrarese and Merritt Ferrarese and Merritt 
(2000) found the black hole mass in bulges to be smaller than that. However our simulations 
start at a situation resembling the final stage of a galactic merger, which means we are 
concentrating on the innermost part of the allready spherical system. 

The black hole particles are initially placed symmetrically about the center of mass of 
the stellar component. Their initial radii are r pa 0.64rh, their initial velocities are 13.6% of 
the circular velocity at this radius. In the given model units, this represents starting points 
for the black holes at x — ±0.5 and v y = ±0.1. The center of mass of the stellar component 
is at the origin. The mass factor between a stellar particle and a black hole particle is: 1338.5 
for 131072 particles, 669.7 for 65536 particles, and 335.4 for 32768 particles. 

In order to have a statistical basis for analysis, we compare the results from five runs 
with 32768 particles, two runs with 65536 particles, and three runs with 131072 particles. 
Not all runs reached the 60 time unit mark due to time step scheduling problems caused 
by accuracy problems in very close encounters between stars and a black hole particle. The 
regularization methods implemented in EuroStar are identical to the ones suitable for open 
or globular cluster simulations. The extreme situation in the late stages of the sinking binary 
black hole problem may cause the chain algorithm to fail (Mikkola and Aarseth 1993). This 
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problem can be solved by applying different regularization methods. However, up to the 
point of failure, the simulations conserved the total energy with relative errors below 10~ 4 . 
In all runs, the binary becomes bound at approximately 10 time units. 

The parameters of the hybrid code have been adjusted in the following way: The SCF 
part uses the parameter a = 0.5 for the basis functions. With this choice, the basis func- 
tions represent a Plummer model to zeroth order, which is in accordance with the models 
used by Clutton-Brock (1973). Also, this choice ensures an optimal representation of the 
actual potential by the expansion method. To allow flexibility in the expansion, seven basis 
functions are used for the radial direction and five (/ = [0...5], m = [—5. ..5]) for the angular 
expansions. The NBODY6++ part uses r\i = 0.01 for the irregular time steps and r\ r = 0.02 
for the regular time steps. Furthermore, the Ahmad-Cohen neighbor scheme (Ahmad and 
Cohen 1973) has been modified in such a way that the search radius for neighbors is en- 
hanced by a factor of 7.7, 14.4, and 27.7 for interactions with the black holes in the runs 
with 32768, 65536, and 131072 particles respectively. 

4.2. The motion of the massive bodies 

In order to compare the runs, all data have been binned by the parameter t, which 
represents the integrated time of the system in A^-body time units (Heggie and Mathieu 
1986). Table 1 gives the number of sample points for the orbital data of the black hole 
particles. For technical reasons, the runs with 131072 particles could not be continued to 
60 time units. When binning the total particle number N tot , table 1 shows the number of 
samples in the row labeled "total". In the following, we present the results for the motion 
of the black hole binary within the stellar system. 

4-2.1. Sinking rate of the binary 

Figure 3 shows the evolution of the quantity (1/a) as the average from the runs above. 
Equation (10) allows us to compute the semi- major axis a from the orbital data of the runs. 
The data is binned for comparison according to the prescription above. From the data in 
each bin we evaluate the average (1/a) and the standard deviation. In order to find the 
hardening rate, we fit a line to the averages, plotted as circles in Figure 3. The standard 
deviations of the data points, given in the plot as the error-bars, supply the weighting factors. 

For the average over all runs plotted in Figure 3, the regression line has a slope of 
8.7 ± 0.4. The dependency on the particle number can be deduced from the data shown 
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Table 1: Number of sample points in the bins used for analyzing the motion of the binary. 
The bins for the evolution time t in iV-body units are centered around < t >bm- Plots using 
bins for the total particle number N tot have the number of sample points given in the row 
labeled "total". 
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Fig. 3. — Evolution of (1/a) as a function of time, after the binary becomes bound. The left 
plot shows the averages computed for each particle group. 32768 data are plotted with open 
triangles, 65536 data with open squares, and 131072 data with circles. The averages over all 
runs carried out are shown in the right plot. The error-bars indicate the standard deviation 
of the data in the bins. 
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on the left side. The slopes are 9.6 ± 0.5 for the 32768 particle simulation, 8.6 ± 0.2 for 
the 65536 particle simulation, and 6.8 ± 0.2 for the 131072 particle simulation. There is 
clearly a dependence of the results for the sinking rate on the particle number. Compared 
with other quantities we analyze in this study, the noise level in the data for 1/a is low. 
We observe strong interactions between the stellar and the black hole particles in runs with 
32768 particles. For this reason, 1/a shows strong steplike changes in both directions. This 
is most likely due to the small particle number. 

4-2.2. Evolution of the eccentricity 

We are only studying the evolution of the eccentricity after the binary became bound. 
Because after 20 time units the eccentricity evolves relatively smoothly for each run, we are 
concentrating our analysis on the time range between 20 and 60 time units. 

Figure 4 shows the mean eccentricity binned in time slots with a width of five time units. 
The symbols represent the averages in these bins. 32768 data are plotted with open triangles, 
65536 data with open squares, and 131072 data with circles. While the eccentricities settle 
at values between e = 0.6 and e = 0.9 for the runs with 32768 particles, the runs with higher 
particle numbers show a fairly parallel evolution. The averages of the 32768 particle runs 
agree very much with the averages from 131072 data. However, the averages for the 65536 
data are clearly higher. 

With our initial conditions, the binary evolves in a highly eccentric orbit, which is 
around e = 0.85. The system retains this high eccentricity until the end of our simulations. 

4-2.3. Evolution of the angular momentum 

In order to study the evolution of the angular momentum of the bound binary we plot 
the angle 9 between the z-coordinate axis and l/l versus time units in Figure 5. 9 is zero 
initially. As with (1/a), all data from the simulations are binned and averaged. The open 
circles in Figure 5 represent the averages, while the error bars are the standard deviations 
in the data. 

Once the binary becomes bound, the 9 changes only slightly in all simulations. Averaged 
over the time between the first bound orbit of the massive particles and the end of the 
simulations, the average value of 9 becomes 0.5 ± 0.3 for the 32768 runs, 0.3 ±0.1 for the 
65536 runs, and 0.5 ± 0.1 for the 131072 runs. 
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Fig. 4. — Evolution of the eccentricity of the massive binary as a function of time. We are 
only plotting the data after the binary has become tightly bound in order to avoid unphysical 
values above 1 and strong scattering of the data. 32768 data are plotted with open triangles, 
65536 data with open squares, and 131072 data with circles. 
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The results in Figure 5 can be fitted by a straight line. The slope of this line is — 0.003 ± 
0.003. When we group the simulations according to particle number, the fitting lines have 
slopes of 0.016 ± 0.004 for the runs with 32768 particles, -0.0044 ± 0.0006 for the 65536 
runs, and 0.006 ± 0.002 for the 131072 runs. Though small, these slopes are all significantly 
nonzero and distinct from one other. Torques clearly act on the binary system throughout 
the simulations. The data for the runs with 32768 particles and with the small mass ratio 
between black holes and stellar particles is very noisy and shows steplike changes in 9. 

While 9 evolves in an ordered way until the binary becomes bound, the angle between 
the rc-coordinate axis and the normalized angular momentum vector behaves more randomly. 
Until the massive particles become bound, <fi changes rapidly reaching all values between 
and 2-7T. However, once the binary becomes bound, settles to a single value for each run. 
The changes in are subsequently of the same order of magnitude as for 9. 

4-2.4- Wandering motion of the binary 

Studying the wandering of the binary using the quantity (r^ om ), we can compare the 
observed motion to the expected Brownian motion in the system. r com is the distance from 
the center of mass of the black hole binary to the center of mass of the stellar system. Figure 
6 implies that the mean motion is not constant with time. However, since the slope of the 
fitting line is (1.0 ±1.1) x 10~ 5 , the behavior is constant within la uncertainty. For the 
individual particle number groups the situation is as follows: For 32768 particles we find a 
slope of (0.6 ± 1.5) x 10~ 5 , for 65536 particles a slope of (1.0 ± 0.7) x 10~ 5 , and for 131072 
particles a slope of (0.9 ±7.8) x 10~ 6 . Compared to its mean value over the whole simulation, 
the evolution of (r^ om ) with time introduces changes of not more than 10%. For this reason, 
we assume (r^ om ) to be constant for the analysis of the Brownian motion. Figure 7 shows 
the mean squared distance between the center of mass of the black hole system and the 
stellar system as a function of the total particle number of the simulations. The slope of the 
fitting line is (—4.5 ± 5.6) x 10~ 9 . Given our small sample of runs we cannot determine a 
dependency of (r^ om ) on the particle number. 

Figures 9 and 8 show the evolution of (v^ om ) as a function of time and total particle 
number N tot . The quantity v com is the velocity of the center of mass of the black holes 
relative to the velocity of the center of mass of the stellar system. For the time dependence 
of (fcom)) we fi n d a slope of (7.6 ± 4.9) x 10~ 5 for the fitting line in Figure 8. For differing 
total particle numbers this slope is (3.5 ± 7.0) x 10~ 5 for 32768 particles, (2.6 ± 2.4) x 10~ 5 
for 65536 particles, and (4.2 ±4.0) x 10~ 5 for 131072 particles. The slope for the dependence 
of (t>com) 011 particle number in Figure 9 is (—1.7 ± 3.9) x 10~ 8 . 
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Fig. 5. — Evolution of polar angle 9 of the angular momentum of the massive binary after 
it has become bound. The angular momentum vector is initially aligned to the z-axis. The 
left plot shows results for the particle groups: 32768 data are plotted with open triangles, 
65536 data with open squares, and 131072 data with circles. The average over all data is 
shown on the right side. 
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Fig. 6. — The mean of r^ om taken over the all simulations as a function of the integrated 
time in iV-body units. The quantity r com is the distance of the center of mass of the black 
hole binary to the center of mass of the stellar particles. The left plot shows results for the 
particle groups: 32768 data are plotted with open triangles, 65536 data with open squares, 
and 131072 data with circles. The average over all data is shown on the right side. 
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number of particles in the simulations. The quantity r com is defined as in Figure 6. 
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4-2.5. Connection between the wandering and the orbital decay 

Figure 10 shows the ratio of the wandering and the semi-major axis of the binary orbit 
a 2 as a function of time. The evolution of this ratio has a strong dependence on the particle 
number, as the wandering is dependent on the simulation size. However, all simulations 
show the same trend in that wandering becomes more important with time for the binary. 
As the right plot in Figure 10 shows, a fitting line with a slope of 1.7 ± 0.6 can fit the data. 
However, the data suggests a nonlinear behavior which should be roughly quadratic, since 
(1/a) increases linearly and (r 2 om ) is roughly constant. 

4-2.6. The effect of dynamical friction 

To study the influence of dynamical friction on the decay of the binary orbit we analyze 
the behavior of its orbital angular momentum as a function of time. As Figure 11 shows, the 
decay shows a two mode evolution. Between and 20 time units, linear regression for (lg(/)) 
gives a slope of (— 6.8±0.9) x 10~ 2 . The line with the more shallow slope (— 1.3±0.2) x 10~ 2 
represents the behavior between 20 and 60 time units. 

4.3. Reaction of the stellar system 

The stellar system reacts to the motion of the black hole in a generic fashion. We find 
that our statistical basis is too small for finding a clear dependency of the results on the 
total number of particles in the simulations. Hence we present only the averages from all of 
our runs. Figures 12 and 13 show the evolution of the density and the velocity dispersion 
respectively for particles within a radius of r csp = 0.032 averaged over all runs as a function 
of time. While the black hole binary becomes bound at ~ 10 time units, the density has a 
maximum at ~ 18 time units, and the velocity dispersion is highest at ~ 23 time units. A 
linear fit (y = a + bx) has been applied to the evolution of p and a between 20 and 60 time 
units as plotted in Figures 12 and 13. For (p), we find a = 2.6 ±0.7 and b = -0.028 ± 0.016, 
for (a) we find a = 2.8 ± 0.6 and b = -0.030 ± 0.013. 
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5. Discussion 

5.1. Hardening rate 

Following Hills (1992), and Quinlan (1996) the hardening rate H of a massive binary 
floating in a sea of light stars is given by, 

^- = H^. (15) 
at a a 

With G = 1 and the assumption that the averages for p and a evolve in the same way 
between 20 and 60 time units, which would render the ratio between p and a constant, 
we find H — 8.7 ± 0.4. This is significantly smaller than the values given by Hills (1992) 
(H = 13.5) and Quinlan (1996) (H « 18). 

Our smaller hardening rate compared to the results of Quinlan (1996) is caused by the 
lower central density and the core type radial density profile of our Plummer model. Quinlan 
(1996) uses Jaffe models for his simulations which allow rapid transfer of orbital energy into 
the dense cusp through tidal interactions. This is also represented in the destruction of the 
cusp Quinlan (1996) observes, while our simulations show a much weaker change for the 
central density. 

Hills (1992) models the shrinking of the binary orbit through three body encounters. 
His greater value of H is consistent with our simulations. We observe steplike changes of the 
binding energy at later times of the simulations, which is less pronounced with increasing 
particle numbers. Because the granularity of the potential is higher in low A^ to t runs, three 
body interactions with the black hole binary become more likely. As shown in Figure 2 such 
three body encounters can enhance the decay of the orbit. Thus, our small H indicates that 
in our simulations shrinking of the black hole orbits is mainly caused by dynamical friction 
and not so much by tidal destruction of cusps or three body encounters. 



5.2. Brownian motion 

If the black holes reach equipartition of kinetic energy with the stars, their expected 
mean square velocity follows from, 

K a ) = ^(vl). (16) 

(Wg qu ) is the mean square velocity we expect for a particle with mass m com , which is the 
combined mass of the two black holes, m* is the mass of the stars and (v%) their mean 
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square velocity. Since our setup involves a Plummer model, we are expecting the binary to 
move in a harmonic potential in later stages of a simulation and for (r 2 om ) oc (v 2 om ). 

While the individual black holes do not reach equipartition, equation (16) can describe 
the Brownian motion of the system. The sum of the black hole masses is represented by 
m com and the center of mass motion of the binary by v com . The mass ratios between the 
individual stars and the black hole binary are 1.49 x 10~ 3 , 7.47 x 10~ 4 , and 3.74 x 10~ 4 in the 
runs with 32768, 65536, and 131072 particles, respectively. However, equation (16) does not 
describe the behavior of the center-of-mass motion correctly. While the velocity dispersion 
drops after 20 time units, the mean square velocity of the black holes (vl) increases. For this 
reason, we compare the measured average (v 2 om ) for the datasets with the mean expectation 
from the right side of equation (16). For (v 2 ) : we take the average over all simulations, which 
is 1.5, neglecting the variability over time. Using this, we can estimate (i>e qu ) and compare 
it with the measured {v 2 om ), 

N tot (v 2 com ) ( 

32768 0.0072 0.0022 3.2 
65536 0.0026 0.0011 2.3 
131072 0.0022 0.0006 3.9 

This means we find a center of mass motion for the binary which exceeds the expected 
value from Brownian motion. However, the motion is enhanced by larger factors than pro- 
posed by Merritt (2000b) or by Chatterjee et al. (2002). A more detailed discussion of this 
result remains for future work. 



5.3. Dynamical friction 

Following Begelman et al. (1980), dynamical friction becomes inefficient as the driving 
forces behind the binary black hole orbit decay after it becomes hard. In order to put 
constraints on this, we estimate the influence of dynamical friction on the decay of the binary 
in the simulation. From Binney and Tremaine (1987), we take the following expression, 
which is derived in Chandrasekhar (1943). 

v. 

, J vU(v*,r)dv* 

= -lQn 2 G 2 log(A) M.(M. + M*) ^ v.. (17) 

at v'l 

Equation (17) describes the deceleration of a particle with mass M, under the influence of 
weak encounters with surrounding particles having a uniform mass M*. In the special case 
of a Plummer model, the integral in equation (17) can be expressed in terms of the escape 
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velocity v esc (Aarseth et al. 1974), 

v. 

l v 2 J(v*, r)dv. = ^ J q 2 (1 - q 2 )Uq, (18) 



where, 

i 

J q 2 (l-q 2 )Uq. (19) 



c 



The quantity n(r) defines the number density of the stars at the position r and q = v/v esc . 
Taking the limit of a continuous system with <C M., the term n(v)M*(M, + M*) becomes 
M.p(r). The integral over q in equation (18) can then be solved in a closed form. 

If the motion of the black holes is determined by their self interaction plus a frictional 
force term, this friction can be linked to the decay of the angular momentum / as follows: 

v 

a = a r + a df — , (20) 
v 

Z = ^(rxv) = ^, (21) 

V V 

where a r is the radial acceleration of the two body motion of the black holes, adf is the 
dynamical friction acting on each black hole, and v is the two-body velocity of the black 
holes. With equations (17) and (18) we can evaluate the impact of dynamical friction on the 
orbital angular momentum according to 



1 16ir 2 G 2 
1 ~ Cvl 



log(A)M.p(r) J q 2 (1 - q 2 )ldq. (22) 



The gravitational constant G is unity in our model units, and the black holes have mass 
M. = 0.01 each. We use the mean orbital velocities for v,. In order to evaluate p(r), we use 
the mean separation of the black holes for r, which introduces only a very small error in a 
Plummer model. Assuming a linear behavior for lg(/) — a + bt, we find /// = &ln(10). Using 
this to estimate the angular momentum from the slopes b of the linear fits in Figure 11, we 
find log(A) ~ 0.15. This result shows that the usual assumption of large A does not hold. 

Both the possibility of a linear fit for the evolution of / and the small hardening rate 
H indicate that mainly dynamical friction causes the shrinking of black hole orbits in our 
simulations. 
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6. Physical units 

As stated before, the collisional simulations which include black hole particles do not 
reach the observed mass contrast in galactic nuclei. In order to transform simulation units 
to physical units, a system size in parsec or a stellar mass in units of solar mass has to be 
chosen. Setting the gravitational constant G, all remaining units can be rescaled (Heggie 
and Mathieu 1986). 

In the following, a run with 65536 particles is scaled to a physical stellar system. Since 
this work focuses on the dynamics of galactic nuclei, the physical mass of the supermassive 
objects motivates the following choices, 

M, = 1.00015 x 1O 7 M , (23) 
M* = 1.49536 x 1O 4 M . (24) 

The masses are chosen so that the total mass of the system is M tot = 10 9 M and the mean 
mass of a particle is M = 1.52588 x 1O 4 M . This choice means that every stellar particle 
with mass M* represents a compact star cluster with the order of 10 4 particles. The chosen 
mass for the black hole particle has approximately the same mass as the central black hole 
of M31 (Magorrian et al. 1998). 

The conversion between physical units and iV-body units follows x p h ys = A conv :r S i m for 
simulated quantities. By choosing the central velocity dispersion to be 110 km/s, we find 

-^conv and -Rconv) 

i? conv = 355.39 pc, (25) 

M conv = lO 9 M , (26) 

T conv = 3.1590 x 10 6 y, (27) 

Konv = HO km/s. (28) 

In a Plummer model, the half mass radius is related to the Plummer radius R by = 
1.30.R (Spitzer 1987). Since R = 37r/16, the half mass radius in model units is r\ = 0.766. 
Therefore, the initial model for the simulated decay of a black hole in the galactic center is 
a Plummer sphere with a half mass radius of 272.23 pc. The initial distance between the 
black holes is 355.39 pc. They become bound after approximately 40 million years. The 
total simulated time is approximately 190 million years. At the end of the simulation the 
black hole distances vary from 1 pc at apocenter to 0.2 pc at pericenter. The semi-major 
axis of the first bound orbit is 21 pc. 

This scaling allows us to compare our results with Begelman et al. (1980). We find 
that our smallest average orbits at the end of the simulation are not yet small enough that 
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gravitational radiation, according to their estimates, would dominate the evolution time 
scale. 

However, at the end of our simulations evolution is still dominated by dynamical friction 
and not by long evolution time scales for hard binaries as proposed by Begelman et al. (1980). 
Their estimate for the gravitational radiation shrinking time scale assumes circular orbits 
for the binary. With eccentricities of roughly 0.85 for the black hole binaries in our runs, we 
expect gravitational radiation to be efficient and coalescence in roughly 10 s years after our 
simulations stopped. 

7. Conclusions 

We have created a new iV-body hybrid code by merging a high accuracy direct Hermite 
integrator of the standard type (Aarseth 1999; Spurzem and Kugel 1999) with a collision- 
less iV-body method which approximates the potential of a given particle distribution by a 
series expansion (Hernquist and Ostriker 1992; Zhao 1996). The SCF method has been 
completely rewritten to include a computation of the time derivative of the gravitational 
force and a fourth order Hermite integrator. We have used this code to model a galactic 
nucleus containing two massive black holes with up to 128k single particles. The evolution 
of the binary black hole is followed from an initial phase, to a phase driven by standard 
dynamical friction where the binary is bound, and then further hardened by three-body 
encounters with single stars. In that hardening phase, we take full advantage of the regu- 
larized three-body integration developed by Mikkola and Aarseth (1996) and Mikkola and 
Aarseth (1998). The method proves to work well, and reproduces standard expectations, 
such as the Chandrasekhar dynamical friction in the initial phase. In the final hardening 
phase due to three-body encounters, we find that the eccentricity of the black hole binary 
maintains a fairly large value (around 0.85). This is very interesting because it decreases 
the time scale for gravitational radiation merger of binary black holes dramatically, thus 
increasing our chances of detecting gravitational radiation from such events with LISA. Due 
to computational limitations, however, our particle numbers are still not large enough to 
fully describe the real physical situation. Any further scaling is problematic, and so further 
work with improved hardware and software must be done. 

We study in detail the motion of a black hole binary in the center of a galaxy. We find 
that the wandering motion does not decay with increasing particle number as expected. The 
mechanism exciting these anomalous motions is unclear. If they exist in simulations with 
realistic particle numbers, they will solve the problem of feeding the black holes with fresh 
stellar dynamical material raised by Gould and Rix (2000). 
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Throughout the computation for the forces and force derivatives in our SCF-scheme 
several special functions have to be tabulated. Recurrence relations provide a very efficient 
means of calculating these functions. The following recursion relations have been applied to 
compute ultraspherical polynomials and their derivatives. As starting values for n 6 0,1 the 
recurrence formulae for the Gegenbauer or ultraspherical polynomials obey the relation 



A. Recurrence relations for Ultraspherical polynomials 




(Al) 



The expressions for higher values are given by: 



2(n + a)^d r r\0-(n + 2a-l)Ci a } l (0 
(n + 1) 



(A2) 



This preprint was prepared with the AAS IATgX macros v5.0. 
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From that the first derivative can be computed as: (Abramowitz and Stegun (1972), equation 
(22.7.22) and table 22.7.) 

r (o+i) f~ _ (n + 2a)ed a) (0-(n+l)Ca(Q ( . 

For practical reasons and higher accuracy the second derivative polynomial is computed 
using equation (A2): 

(Q+2) _ 2 (n + a + 2)£Ct +2) (Q - (n + 2(q + 2) - 1) d«f } (Q ... 

G " +1 (^Ti) ( j 

Because the particle track is approximated by using the Hermite scheme, one has to find 
forces and the first force derivative simultaneously. An approximation using two timesteps 
for the first force derivative introduces errors to the second and third derivative of the forces. 
All particles move within a time dependent potential; therefore, the first derivative has a term 
describing the change of the potential and a term describing the change of force depending 
on the particle's orbit. 

With the help of the orbit integration for the single particle in a given static potential case, 
equation (A5) evaluates to: 

d da r da r da r ■ da r . ■ . . 

- a(t, r) = (— + — r + — * + — <p - atf - w sm tf) e r 

da, & da# dan ■ da# ■ . 



da^ da v da v k da v . . . , 

+ (^rr + o r + ~^T^ + o V 9 + a r <p sm-d + a$(p cos ■&) e v 
at or ov d(p 



The evalutation of the first term on the right hand side of equation (A5) is given in 
section B. The derivatives with respect to the spatial coordinates in equation (A6) can be 
found in section C. 



B. Time-dependency of the potential 



Because all positions and velocities of the dataset are time- dependent, the partial deriva- 
tives with respect to t apply only to the coefficients A n i m . These are implemented as the vari- 
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ables Ci m (r), Di m (r), Ei m (r), and Fi m (r), from which the partial derivative can be formed: 
g l 7 = N lm A nl ® nl {r) m kQ- t {®m{r k ) P lm {cos{-d k )) cos(m^ fe )J , (Bl) 

n=0 k 



3D 



= JV, m ^ AA(r) mk Qt (®ni(r k ) P lm (cos(# k )) sin(m^ fc )J , (B2 

n=0 fc 

= N lm jr A nl J2 m "§i (® nl{rk) Pl ^ cos ^ <x*(™Pk)) , (B3) 

n=0 fc 

= N lm Y,A nl ^—J2rn kFt {* 

n=0 k 



dt 



'ni(r k ) Pi m (cos(i3 k )) cos(my? fe )) , 
>ni{r k ) Pi m (cos sin(my9 fe )) . (B4) 



With: 



(9 ~ ~ 



j rj 2/ + 1 | 4r° a(2Z + l) + i C^fo) 



— P Zm cos [d k )) = — — sin tf fc — o 7^ — , 

OT ot acos(i?fc) 

d 



cos(mip k ) = - m -T^- sm(mip k ), 
^ sin(my9 fc ) = rn^- cos(mip k ). 



(B5) 

(B6) 

(B7) 
(B8) 



in the coefficient computation section the standard leap frog integrator provided by 
Hernquist and Ostriker (1992) is extended by two additional variables, which are computed 
by using the recursion relations in section A. 



C. Orbit dependency of the force derivative 



In order to account for the change of force due to the particle orbit one has to calculate 
the nine partial derivatives in equation (A6). These nine derivatives will now be listed. In 



order to save some space the second derivative of $ n z( r ) is given first: 



d 2 ~ 

— <$> nl {r) = $„,(r) 



I r~ (2Z + 1) 
r 

(21 + 1) r 



I 



r r (l + r « 
i_ 

a — 1 + ar a ) 



(l + ri)2 ar 2 



a; 



r a(l + r-' 



+ 



- — (f — a — (a + l)r«) 



/ ro (2/ + 1) 
r r (l + r«) 



r 2 a 2 (l + r « ) 3 

Vary' (l + r ^)4 #(0 
The nine derivatives can be implemented as follows: 



#(0 



C.l. Derivatives with respect to r 

The radial derivative for the radial acceleration becomes: 
da r 

= -2 

=0 m=0 



dr 



- ^2 S ^™( cos (^)) [ G «m( r ) cosfmy?) + if im (r) sin(my9)] , 



with: 



d 2 „ 



G/ m = N lm ^ A n i—$ n i(r) ^ m k $ n i(r k ) P lm (cos(-& k )) cos(rrup k ), 



n=0 

oo 



' dr 2 

d 2 ; 



= ^2 Anl ~Q^2® nl ( r "> ^2 m k^nl(r k ) Plm(cOs(l3 k )) COs(my9 fc ) 

n=0 k 

The radial deriavtive for the acceleration in d direction becomes: 

oo oo 



da# ■ /cnV^ <9P im (cos(^)) 

= - sin i!' > > — — — 

^ ^ <9cos($) 

1=0 m=0 v 1 



dr 



x 



\c lm (r) - -E lm (r) ) cos(my?) 



lAmW - -Fi m (r) ) sin(my?) 
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The radial deriavtive for the acceleration in cp direction becomes: 



dr 



mPi m (cos(fl)) 
sin(tf) 



oo oo 

EE 

1=0 m=0 

( \D lm (r) - -F lm (r) ) cos(/;v) 



^Cim(r) - -E lm (r) ) sin(my?) 



(C6) 



C.2. Derivatives with respect to d 

The derivative with respect to ■& for the radial acceleration becomes: 



da r 



sin 



w E E ^l c ff )} [E ' m(r) cos(mv?) + Flm{r) sin(mv?)] 



2=0 rrt=0 

The derivative with respect to $ for the acceleration in i? direction becomes: 



da<> 1 ^ ^ / <9P m (cos($)) . 2 d 2 Pi m {cos(d)) 

W = r E E g cos( ^) - Sin W gcosW^ 

2=0 m=0 v v 7 v 7 

x [C; m (r) cos(mv9) + £) im (r) sm(m</?)]. 
The derivative with respect to $ for the acceleration in <p direction becomes: 
da v f dP lm {cos{$)) P im (cos(tf))\ 

w = r E E - { dcoaW + cos W sin2w J 

2=0 m=0 \ > / 

x [D ;m (r) cos(mv9) - Cz m (r) sin(m</?)]. 
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(C8) 
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C.3. Derivatives with respect to </? 

The derivative with respect to tp for the acceleration in radial direction becomes: 

q OO OO 

Tp = - ^2 E m ^™( C0S W) [^m(^) cos(my?) - P im (r) sin(m<p)]. (CIO) 

^ 2=0 m=0 

The derivative with respect to <p for the acceleration in i? direction becomes: 

da* sin(tf) ^ ^ dP; m (cos(tf)) . 

2^ m — — [AmW cosfmy?) - C lrn (r) sm(mip)\. (Cll) 



dcp r ^ 

T 2=0 



rn=0 



<9cos(-#) 
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The derivative with respect to ip for the acceleration in p direction becomes: 

^ = \ E E ^ [<W0 cos (^) + A m (r) sin(m^)]. (C12) 

^ Z=0 m=0 ' ' 
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Fig. 8. — The mean of v^ om as a function of the integrated time in iV-body units. The 
left plot shows results for the differen particle groups: 32768 data are plotted with open 
triangles, 65536 data with open squares, and 131072 data with circles. In the right plot 
the data for v^ ora have been averaged over all particle groups, the error-bars represent the 
standard deviation in the data. 
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Fig. 9. — The mean of v% om as a function of the total particle numbers in the simulations. The 
values for v^ m have been averaged over the total simulation time, the error-bars represent 
the standard deviation in the data. e eom is the relative motion of the center of mass of the 
massive binary relative to the center of mass of the stellar system. 
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Fig. 10. — Wandering of the binary in relation to the squared semi-major axis a 2 of the bound 
black hole binary as a function of time. The left plot shows the results for each particle group: 
32768 data are plotted with open triangles, 65536 data with open squares, and 131072 data 
with circles. The right plot shows the results for all simulations put together. 
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Fig. 11. — The evolution of the orbital angular momentum as a function of time for the 
collected data of the runs. The left plot shows the averages computed for each particle 
group. 32768 data are plotted with open triangles, 65536 data with open squares, and 
131072 data with circles. The right plot shows the averages for all particle groups together. 
The error bars represent the standard deviation in the data. In order to distinguish between 
the two modes of evolution, linear regression is applied to the bins between and 20 and 20 
and 60 time units separately. 
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Fig. 12. — Evolution of the stellar density in a central sphere of the cluster with r csp = 0.032. 
p is the average over all simulations, the error-bars represent the standard deviation in the 
data. The dot-dashed line shows our linear fit for the evolution between 20 and 60 time 
units. 
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Fig. 13. — Evolution of the stellar velocity dispersion in a central sphere of the cluster with 
^cs P = 0.032. The quantity a is the average from all simulations, the error-bars show the 
standard deviation in the data. The dot-dashed line shows our linear fit for the evolution 
between 20 and 60 time units. 



